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PACS 42 . 25 . Dd - Wave propagation in random media 

PACS 42 . 68 . Ay - Propagation, transmission, attenuation, and radiative transfer 

PACS 05. 40. -a - Fluctuation phenomena, random processes, noise, and Brownian motion 

Abstract. - Two-frequency radiative transfer (2f-RT) theory is developed for geometrical optics 
in random media. The space-frequency correlation is described by the two-frequency Wigner 
distribution (2f-WD) which satisfies a closed form equation, the two-frequency Wigner-Moyal 
equation. In the RT regime it is proved rigorously that 2f-WD satisfies a Fokker-Planck-likc 
equation with complex-valued coefficients. By dimensional analysis 2f-RT equation yields the 
scaling behavior of three physical parameters: the spatial spread, the coherence length and the 
coherence bandwidth. The sub-transport-mean-free-path behavior is obtained in a closed form by 
analytically solving a paraxial 2f-RT equation. 
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• i-H Introduction. — Correlation functions of fields arise 
^naturally in the description of fluctuations and are ubiq- 
^^uitous objects in statistical physics. The most basic of 
Cl^hose are the second-order correlations in the space-time 
or space-frequency domain; the two are equivalent to each 
other via the Fourier transform. When the field fluctua- 
^►tions can be described as a Gaussian stochastic process, 
OSall the correlation functions of the field can then be ex- 
0^ pressed in term of the second order ones, by the use of 



the moment theorem for Gaussian processes. The second 
order space-frequency correlation then emerges as an in- 
dispensable tool for studying fluctuations of fields and is 
equivalent to the mutual coherence function describing the 
infield correlation at two space-time points [1]. 

^ Spatial and temporal structures of ultrawide-band high- 
• ^ frequency fields can be appreciably affected by small ran- 
rSdom changes of the medium parameters characteristic of 
^ almost all astro- and geophysical environments. An impor- 
tant step toward analytical understanding of pulse prop- 
agation in multiply scattering media is then to derive the 
equation for the space-frequency correlation, obtain the 
qualitative information about its behavior and, if possible, 
find its (asymptotic) solutions. This problem has been ex- 
tensively studied in the literature, see, e.g., [3-7]. The 
main distinction of our approach from previous works is 
that our approach to space-frequency correlation is carried 
out in terms of the two-frequency Wigner distribution (2f- 
WD) for which we will derive rigorously equations of rel- 



atively simple form in the radiative transfer (RT) regime 
and obtain an exact solution for the small-scale behavior 
below the transport mean-free-path [1,8]. 

The standard (equal-time or -frequency) Wigner dis- 
tribution (WD) is a quasi-probability density function in 
phase space and was first introduced by E. Wigner [9] in 
connection to quantum thermodynamics and later found 
wide-ranging applications in classical [10], [11], as well as 
in quantum optics [1], [12]. In classical optics, a main use 
of the Wigner distribution is connected to high-frequency 
asymptotic and radiative transfer, both of which can be 
most naturally worked out from the first principle in phase 
space (see the review [13], [14] and references therein). 

The main advantage of 2f-RT over the traditional equal- 
time radiative transfer theory is that it describes not just 
the energetic transport but also the two space-time point 
mutual coherence in the following way. 

Let the scalar wave field Uj,j = 1,2, of wavenumber 



j • ■ 



1,2. be governed by the reduced wave equation 



AU j (r) + k*(v + V(r))U j (r) = 0, rel 3 , j = l,2 (1) 

where v and V are respectively the mean and fluctuation 
of the refractive index associated are assumed to be real- 
valued, corresponding to a lossless medium. For simplicity, 
we restrict our attention to dispersionless media (see [2] 
for discussion on dispersive media). Here and below the 
background wave speed is set to be unity so that kj — ujj . 
Let u(tj,Xj),j = 1, 2 be the time-dependent wave field at 
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two space-time points (tj,Xj),j = 1,2. Let x = (u)\Xi + 
w 2 X2)/2 and y = W1X1 — cj 2 x 2 . Then we have 



(u(ti,Xi)u*(t 2) X2)) (2) 
J e ip-y e i(«2ta-a,iti) (W(x,p;o;i,W2)>dwidw2dp 



area, A the wavelength of radiation and L the distance 
of propagation. Let Nf = XL/A be the Fresnel number 
and the transport mean-free path. The dimcnsionlcss 
conductance can then be expressed simply as g = k£*/Nf. 
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where W(x, p;wi,w 2 ) is the 2f-WD defined by 



W(x,p;wi,W2) 
1 

(27T) 



With the scaling (28), M 
9 ~ Q- 2 £- 2 > 1 for 0e < 1. 

To adapt to the weak coupling and the geometrical op- 
tics (see below) scalings we introduce the two parameters 
e, 9 into the 2f-WD and redefine it as 



- tp - y Ui(- + ^)u;(--^)dy 



W s (x,p) 



(4) 
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(2tt)3 ./ ' ""/,-, 2k 1 )U * i k 2 2k 2 )dy 



Here and below (•) denotes the ensemble average. For 
temporally stationary signals, wave fields of different fre- 
quencies are uncorrelated and only the equal-frequency 
WD is necessary to describe the two-time correlation. In 
comparison, the single-time correlations with t\ = t 2 = t 
gives rise to the expression 

(u(t,xi)u*(t,x 2 )) = J div'dp e tp - y e luJ,t 

x J du){W{x,p;u- w'/2,u + u'/2)} 

which is equivalent to the central-frcqucncy-intcgrated 2f- 
WD. For ease of notation, we will drop the frequency ar- 
guments when writing the 2f-WD below. 

Weak coupling limit. — The radiative transfer 
regime sets in when the scale of medium fluctuation is 
much smaller than the propagation distance but is com- 
parable or much larger than the wavelength. Based on the 
general principle of central limit theorem, RT corresponds 
to the scaling limit which replaces v + V in eq. (1) with 



^(u+VeV( r -)), 9>0, e«l 



(3) 



where e > and _1 > are, respectively, the ratio of 
the scale of medium fluctuation to the propagation dis- 
tance and the wavelength. Thus, s9 is the ratio of the 
wavelength to the propagation distance and as a result 
we rescale the wavenumber as k — > k/(s9), giving rise to 
the prefactor (9e)~ 2 . This is so called the weak coupling 
(or disorder) limit in kinetic theory [15] under which the 
Anderson localization can not take place. 

We assume that V(x) is an ergodic, mean-zero, statis- 
tically homogeneous random field. As a consequence, V 
admits the spectral representation V(x) — J e tx ' p V(dp) 

where the spectral measure V satisfies (v(dp)V(dq)^ = 

S(p + q)<j>(p)dpdq with $ is the power spectral density. 
Since V is real-valued, <&(p) is real-valued, non-negative 
and possesses the symmetry ^(p) = <]>(— p),Vp. 

Physically speaking radiative transfer belongs to the dif- 
fusive wave regime under the condition of a large dimcn- 
sionlcss conductance g » 1. Let A be the illuminated 



In view of the definition, we see that both x and p are di- 
mcnsionlcss. The particular scaling factors are introduced 
in (4) so that W £ satisfies the following Wigncr-Moyal 
equation exactly [2] 



p • VW E = -^LW e 
where the operator C is defined by 

cw^{x,p) = -Ljv { d q y^wc(x,p-^-) 

9q 



(5) 
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J v(d q y^w%x, P + ^). 



In contrast, the Sudarshan equations for the mutual coher- 
ence function are, like (5), first-order in time but nonlocal 
in space even in the case of free field [1]. 

High-frequency regime. — Before we consider the 
radiative transfer limit e [ further let us take the high- 
frequency limit 9 { while maintaining the following rela- 
tionships 



lim ki — lim fc 2 

0^0 0->O 

fc 2 - fci 



9ek 



= k 
= 



(6) 



where (3 > is independent of 9 and s, representing the 
normalized difference in wavenumber. Frequencies within 
the range described by (6) remain coherent with one an- 
other. 

In this regime, we see from (4) that to leading order 
the center of two field points is x/k and the difference is 
9e(y + (3x)/k. Passing to the limit 9 1 in (5) we obtain 
the first-order partial differential equation 



p-V x W £ (x,p) 

1 fx 

~2ivf( w H^ 



(7) 



V p -ipx W £ {x,p). 



For (3 = 0, eq. (7) is the static Liouville equation. For 
/3 > 0, eq. (7) retains the wave character and is the focus 
of the subsequent analysis. We shall refer to eq. (7) as the 
two- frequency Liouville equation (2f-LE). 
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Consider, for instance, the WKB ansatz Convergence to diffusion in momentum. — If V 

decorrelates sufficiently rapidly (see [16] for a precise for- 

Uj(r) = Aj(r) exp (^Sj(r)\ , j = 1, 2 mulation), the probability distribution of (x £ (-), p e (-)) de- 

^ " e ' fined by (14)-(15), converges weakly, as e — > 0, to that of 

where the phase Sj and the amplitude Aj depend on the thc Markov process (x(-),p(-)) where 

frequency differentiably. In the first case, assume S\ = r-t 

S 2 = S. Then in the high-frequency limit 2f-WD becomes x(i) = x - / p(s)ds (16) 

Jo 

IE £ (x,p) = e i0x ' p e~ i0ks{ ^\A\ 2 (^j8(p - VS(-)) (8) and p(-),p(0) = p, is a diffusion process generated by thc 

' operator 

which satisfies 2f-LE. In the second case, assume Sj(r) = \ 

kj ■ r, |kj| = 1, with thc additional conditions — 4fc^ p ' ^ ' ^ p 

lim k lim k k (9) w ^ ^ e ( momen t um ) diffusion coefficient 

k 2 -ki A1 « D(p) =tt / $(q)<5(p-q)q®qdq. (17) 

- Ak (10) 7 

Writing D as 

where Ak is independent of 8,e. Then the the 2f-WD r 

becomes D(p) = tt / $(q)<5(p • q)II(p)q ® II(p)qdq (18) 

|A| 2 ^ — ^ e^* 4 ' x (5(p — k) (11) where II(p) is the orthogonal projection onto thc hyper- 
plane perpendicular to p we see that the momentum dif- 

where (3 is absent due to the linear phase profile Sj. fusion process is concentrated on the sphere of radius |p|. 

Given, say, the Dirichlet boundary condition F imposed In other words, the limiting kinetic energy |p(i)| 2 /2 is pre- 

on the boundary dV of a phase-space domain T>, 2f-LE can served by the elastic scattering process. This observation 

be solved by the method of characteristics as shown below, will be useful for the subsequent calculation. 

Thc form of 2f-LE suggests the "gauge transformation" of The consequence is the convergence of the ensemble av- 

2f-WD erage (2» £ (x,p)) to 

2U e (x,p) = e- 4 ' 3x - p V^ £ (x,p) (12) 2U(x,p) (19) 

= E ( e - i ' 3 lpl 2T - i / 3x ( r )-p( T )^(x(T) p(t))1 

which then satisfies the following more convenient x p I v v ' v "} 

equaiton where r = r(x, p) is the first passage time of the Markov 

p • V X 2U £ + i(i\p\ 2 W e (13) P rocess ( X W> P(*)) with x (°) = x > P(°) = P and E x,p thc 

x 11 corresponding average. 

= (^^) ( — ) • V P 2U £ Now let W(x, p) be the solution of the following bound- 

2k^/e \ekl ary yalue prorj i em: 

with the boundary condition that 2U £ (x,p) = 1 / \ / , fl \„, 

i -a lev < vt \ t t P • V x ff = — V p — i/3x ) • D • V p — ifix I W (20) 

exp [— i(3x ■ p]F(x,p) = tf(x,p) on oD. In view of 4fe V / \ v ) 

(12) 2B £ is the Fourier transform of the two-point with w = p on Qv and wc wU1 ghow that thfi ^ 

function L7i <8> t/ 2 * in the location difference (i.e. y + /3x mtion of 2f _ RT ig the pointwisc limit of the avera ge 2f- 

measured in the unit of the central wavelength). WD Eq ^ (2Q) ig our two-frequency radiative transfer (2f- 

Consider the Hamiltonian system of time-reversed char- RT) equation Because we have considered thc high frc- 

acteristic curves quency asymptotics the scattering term takes thc form of 

d a second-order differential operator rather than the more 

—x (t) — — p (t) (14) f arm ]j ar integral operator. 

^ I /x £ (t)\ Let p(t) be the diffusion process generated by the gen- 



Jt P%t) = -_(VV)(— j (15) era tor A and define 

W(i,x,p)=cxp -i(3t\p\' 2 -if3x-p W(x,p). (21) 



with x £ (0) = x,p £ (0) = p. Let t £ = r £ (x,p) be the first 
passage time when the trajectory (x £ (-), p £ (-)) hits thc By Dynkin > s f ormu i a [17] we have that 
boundary of the phase-space domain T>. We then have 



E 



{lE(r,x(r),p(T))} = TE(0,x,p) 



2» £ (x,p) 

= e -^/ T£ |p £ ( S )| 2 d S -^(r=).p^(^) F(x£(T e ))p£(r£)) _ +Ex 'p{/ [^-p-V x + ^]W(s,x(s),p(a))d«}. 



p-3 



A. C. Fannjiang 



From (20)-(21) it follows that 



d 



[^-P-V x + .4]^ = 



and 



o~ x <J P ~ 1//3 C whose left hand side is the space-spread- 
bandwidth product. Finally the removal of the constant 
C determines a p from which a x and f3 c can be determined 
by using the duality relations. We obtain 

a p ~ fc 2 / 3 ^ 3 , a x ~ k-^C 1 / 3 , (3 C ~ fc 2 / 3 C- 2 / 3 . (28) 

2U(x,p) = E x , p {t^(t,x(t),p(t))} = e- l/3xp VF(x,p). (22) Spatially anisotropic media. - Forward-scattering 

approximation, also called paraxial approximation, is valid 
when back-scattering is negligible and, as we show below, 
this is the case for spatially anisotropic media fluctuating 
slowly in the (longitudinal) direction of propagation. 

Let z denote the longitudinal coordinate and xj_ the 
transverse coordinates. Let p and p^ denote the longitu- 
dinal and transverse components of p £ K 3 , respectively. 
Let q = (q, qj_) £ R 3 be likewise defined. Consider now 
a highly anisotropic spectral density for a medium fluc- 
F(x(r), p(r)) 1 (24) tuating much more slowly in the longitudinal direction, 
i.e. replacing ^(q) in (17) by 7y _1 <I> (r)~ 1 q, qj_) , rj <C 1, 
which, in the limit r\ — > 0, tends to 



Therefore, in view of (12), W(x, p) is the pointwise limit 
of (W e (x, p)). It is straightforward to check that 2U is the 
solution to the equation 



p V x W+if3\p\ 2 W = AW. 



(23) 



From (19) and (22) we obtain the probabilistic repre- 
sentation for W 



^(x,p) = E x , p { e - 



i/3/ T x( s )dp( s ) 



where J* x(s)dp(s) is understood as an Ito integral [17]. 
Expression (24) suggests a numerical solution procedure 
for 2f-RT by Monte Carlo simulation. 

Isotropic medium. — Eq. (23) clearly is translation- 
ally invariant in x due to the stationarity of the medium. 
If the medium is also statistically isotropic, then eq. (23) 
is rotationally invariant. To see this let us consider an 
isotropic spectral density ^(p) = <j>(|p|). Then we have 
D = C|p|- 1 n(p) where 



°= I J <W Si >*«<>""> 



(25) 



is a constant. The coefficient C (and D) has the dimension 
of inverse length while the variables x and p are dimen- 
sionless. 

The resulting A is invariant with respect to rotation in 
p. Hence if 2U(x, p) is a solution to (23) then 2U(i?x, Rp) 
is also a solution where R is any orthogonal matrix. 

Spatial (frequency) spread and coherence band- 
width. — Through dimensional analysis, eq. (20) yields 
qualitative information about important physical param- 
eters of the disordered medium. For this, let us assume 
for simplicity the isotropy of the medium as above. 

Now consider the following change of variables 



x = cr x fcx, p = a p p/k, (5 = f3 c j3 



(26) 



where a x and a p are respectively the position spread and 
the spatial frequency spread, and (3 C is the coherence band- 
width, also known as the Thouless frequency. Let us sub- 
stitute (26) into eq. (20) and aim for the normalized form 

p . V X IF = (Vp - ifa) ■ M (V p - i0Z) W. (27) 

The 1-st term on the left side yields the first duality re- 
lation <T x /a p ~ 1/fc 2 . The balance of the terms in each 
pair of the parentheses yields the second duality relation 



S(q) J dwQ (to, qj_). 



(29) 



We then obtain the transverse diffusion coefficient 



Di(pi) 



dw$(w, q±)5(p± ■ q±)q± ® q± 



whereas the longitudinal diffusion coefficient now vanishes. 
In other words, the longitudinal momentum is decoupled 
from the transverse momentum and is not directly affected 
by the multiple scattering process. 

For simplicity we assume the transverse isotropy, i.e. 
$(w,P_l) = $(w,|p_l|), so that D ± = C^p^U^p^) 
where 



C ± 



|P± I 



|qx| 



®(w, \<l±\)h±\dwdq A 



is a constant and IIj_(pj_) is the orthogonal projection 
onto the line perpendicular to pj_. Hence eq. (20) reduces 
to 



pd z + p ± ■ V Xi Jl4> 

Cj - m -a \ n x(Px) m 
- ^(Vp,-,/3x ± ).^^(V p , 



(30) 



i/3x_l) W. 



Note that the longitudinal momentum p plays the role 
of a parameter in eq. (30) which then can be solved in 
the direction of increasing z as an evolution equation with 
initial data given at a fixed z. 

As before we can obtain the scaling behaviors of spa- 
tial spread, coherence length and coherence bandwidth by 
dimensional analysis. Let cr* the spatial spread in the 
transverse coordinates x^, l c the coherence length in the 
transverse dimensions, f3 c the coherence bandwidth and L 
the distance of propagation. We then seek the following 
change of variables 







xx = -^r, p ± = p ± M c , *=4r, , 



(31) 
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to remove all the physical parameters from (30) Follow- 
ing the same line of reasoning, we obtain that £ c er* ~ 
L/k, cr*/4 - 1//3 C , 4 - C~ 1/3 L-V3fc-i and hence 



Ci /3 L 4 / 3 , 



C ~ crr^L-Wk- 1 . 



Small scale asymptotic. — On the scale below the 
transport mean-free-path the scattering is extremely 
anisotropic and the scattering amplitude is highly peaked 
in the forward direction. This observation leads to a 
paraxial approximation of 2f-RT which turns out to be 
analytically solvable. 

Let z be the direction of propagation of a collimated 
beam. On the scale below the 2f-WD near the source 
point would be highly concentrated at the longitudinal 
momentum, say, p = 1. Hence we may assume that the 
projection Il(p) in (18) is effectively just the projection 
onto the transverse plane coordinated by Xj_ and we can 
approximate cq. (20) by 



where 



- J <&(0,q ± )|qj 2 dq ± . 



Here we have assumed, for simplicity, the transverse 
isotropy of <i>. Eq. (32) is another form of paraxial approx- 
imation for which only the one-sided (incoming) boundary 
condition (at z = const.) is needed. 

We use the change of variables (31) with l c ~ 

k ^ Gj_ ^ ' 00 ^* G A_ ' f^ c ^ ^ ^* ^ tg- 

move the physical parameters from eq. (32). The trans- 
port mean-free-path can be determined by setting £ c ~ 
l,i.e. £* ~ fe- 2 Cj x . Perfor ming the inverse Fourier trans- 
form in p on the rescaled equation we obtain 

d s r-iVy ± Vx x r = -| y± + /3x ± | 2 r (33) 

which is the governing equation for the two-frequency co- 
herence r. By a simple change of coordinates, eq. (33) 
can be converted into a form similar to the time dependent 
Schrodinger equation with a (purely imaginary) quadratic 
potential and then solved analytically. Let Ar = yj_+/3xj_ 
and Ar' = y' ± +f3x.' ± be the field point offset and the source 
point offset, respectively, measured in the unit of central 
wavelength. The propagator for the initial value problem 
from the source point (xj_, Ar) to the field point (x^, Ar') 
is given by [2] 



(14/3) 1 / 2 



(27r) 2 zsinh [(i4/3)V2£] 

Ar- 



_ e Ar ~ ~ Ar ' 



(i4^)!/2 



ah (i4^) 1 / 2 £ 



xe 



(•id) 



172- 



Ar' 




5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 



Fig. 1: The absolute value of (34) as a function of z £ [0.5, 1] 
for Ar = Ar' = 1,/? = 0.3, 1,3.3 in solid, dashed and dotted 
lines, respectively. 

which converges, in the limit (3 [ 0, to the propagator for 

(3 = 



i^fW (32) { 2wz)~ 



2 e f(x x -xl)-(Ar-Ar') e - 



KIAi-r+Ar-Ar'+IArY 



(35) 



The quadratic-in-Ar nature of the exponents appearing 
in (34)-(35) is the consequence of the paraxial approxima- 
tion. Expression (35) is related to the asymptotic solution 
of the Schwarzschild-Milne equation in the case of very 
anisotropic scattering [18]. 

In view of (9)-(ll), to get the the correlation of two in- 
cident plane waves we simply express (34) in the variables 
x_l,x^_ and yj_,y^_ and integrate it with e lAk x i e jk y-L. 

The functional form of (34) in its dependence on (3 and 
z is the main characteristic of the sub-f*-scale behavior 
(see Fig. 1). 

Conclusion and discussion. — The main contribu- 
tion of the present Letter is the rigorous derivation of the 
2f-RT equation (20) governing 2f-WD in disordered media 
and the probabilistic representation (24). As a result, by 
(2) we can express the two-space-time correlation as 



(u(£i,xi)u*(£ 2 ,x 2 )) 



Jpy Jk(t 2 -t 1 )Ak6el3(t 1 +t2)/2 



W(x,p)dkd/3dp 



(34) 



where 20 is the solution to eq. (23). The medium 
characteristic enters the Fokker-Planck-like eq. (20) only 
through the momentum diffusion coefficient (17). By di- 
mensional analysis with (20) and its variants we obtain 
scaling behavior of spatial spread, coherence length and 
coherence bandwidth for isotropic and anisotropic me- 
dia. We also show that the paraxial regime is valid for 
anisotropic scattering, giving rise to two forms of parax- 
ial 2f-RT equations. Finally by solving one of the parax- 
ial equation (32) we obtain precise profile of the space- 
frequency correlation on the scale below the transport 
mean-free-path. 

Let us compare our results, especially (34), with the 
existing results in the literature which mostly concern with 
the bulk behavior of the space-frequency correlations. 
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Since the bulk behavior concerns the scales larger 
than the transport mean-free-path the existing results are 
mostly based on the diffusion approximation to the dis- 
placement process x(i) or the random-matrix method (see. 
e.g., [5], [19] and references therein). The diffusion regime 
represents an isotropic scattering under the condition of 
cquipartition of energy while the small-scale asymptotic 
(34) describes an extremely anisotropic scattering. 

Clearly the diffusion approximation is unsuitable for 
evaluating (24) because of the presence of the ltd integral 
with respect to the momentum process p(t). Therefore to 
get the two-frequency coherence, the notion of the interfer- 
ence of diffusions is invoked via diagrammatic techniques, 
see the review [5]. 

In the diffusion approximation for isotropic media, the 
(dimentionless) x-diffusion coefficient can be derived 
from (20) with (3 = 



4Arjp>| 
3C 



(36) 



The momentum-dependence of JJ* shows that even in 
the diffusion approximation the momentum variable is 
only hidden from sight. With (36) and (28) we can 
rewrite the scaling behaviors of the spatial spread, the 
spatial frequency spread and the coherence bandwidth as 
o x ~ k~ 1 D« 1 ^ 3 , a p ~ kD* 1 ^ 3 ,{3 C <~ D^ 3 . 

The short-range correlation C\ of wave intensities 
propagating through disordered media is manifest in 
the speckle pattern. C\ can be obtained by squar- 
ing the two-frequency coherence of the wave fields [20] 



and the commonly accepted form is exp [— 2y2/3] which 

is just the large (3 asymptotic of the squared factor 
| sinh [(M/?) 1 / 2 ^]!- 2 at z = 1 (see, e.g., [21-23]). 

More precisely, the squared absolute value of (34) for 
z = 1 and median to large (3 is approximately given by 



(27T) 



|Ar|" 



/2f> g 



(37) 



Expression (37) is essentially the same as the parax- 
ial approximation of the short-range correlation C\ re- 
viewed in [5]. The multiplicative nature of (37) 's func- 
tional form in Ar and Ar' is consistent with the same 
structure in the short range intensity correlation C\ = 
A(Afc)F(Ar)F(Ar') discovered in [24]. Again, the Gaus- 
sian form in (37) is different from the form-factor F in [24] 
due to the paraxial approximation made in obtaining (37) . 

The long- and infinite-range correlations, represented 
by C2 and C 3 respectively, can also be obtained by our 
method, [5,24-27]. The calculation is much more involved 
and will be presented elsewhere. 
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